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Abstract 

A translation invariant A^-polaron system is investigated at arbitrary electron-phonon coupling 
strength, using a variational principle for path integrals for identical particles. An upper bound for 
the ground state energy is found as a function of the number of spin up and spin down polarons, 
taking the electron-electron interaction and the Fermi statistics into account. The resulting addition 
energies and the criteria for multipolaron formation are discussed. 
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I. INTRODUCTION 



Thermodynamic and optical properties of interacting many-polaron systems are intensely 
investigated, because they might play an important role in physical phenomena in high- Tc 
superconductors (see, e.g., Refs. [l|, 0] and references therein). In particular, numerous 
experiments on the infrared optical absorption of high-T^ materials (see, e.g., Refs. 0,0,05 
EIqII]) reveal features which are associated with large polarons 0,11,0]. 

For the case of weak electron-phonon coupling strength, a suitable variational approxima- 
tion to the ground state energy of an interacting many-polaron gas was already developed in 
jlo| , using a many-body canonical transformation for fermions in interaction with a phonon 
field. The static structure factor of the electron gas is the key ingredient of this theory. 
Based on the approach of Ref . [3| , a many-body theory for the optical absorption at a gas 
of interacting polarons was developed jll|. The resulting optical conductivity turns out to be 
in fair agreement with the experimental band" by Lupi et al. in the optical-absorption 
spectra of cuprates. 

At arbitrary electron-phonon coupling strength, the many-body problem (including 
electron-electron interaction and Fermi statistics) in the A^-polaron theory is not well de- 
veloped. Within the random-phase approximation, the optical absorption of an interacting 
polaron gas was studied in Ref. ^^l? taking over the variational parameters of Feynman's 
polaron model which however are derived for a single polaron without many-body ef- 
fects. For a dilute arbitrary-coupling polaron gas, the equihbrium properties [ij, [l^ and 
the optical response jl^ have been investigated using the path-integral approach taking into 
account the electron-electron interaction but neglecting the Fermi statistics. Recently, the 



formation of many-polaron clusters was investigated in Ref. [17| using the Vlasov kinetic 
equations However, also this approach does not take into account the Fermi statistics 
of electrons, and therefore it is o nly valid for sufficiently high temperatures. 

The path integral treatment [i^, 21 1 of the quantum statistics of indistinguishable 



particles (bosons or fermions) provides a sound basis for including the many-body effects in 
a system of interacting polarons i^l . This approach was used |23l . 0, 2^ for calculating 
the ground state energy and the optical conductivity spectra at arbitrary electron-phonon 
coupling strength for a finite number of interacting polarons in a parabolic confinement 
potential. However, the translation invariant polaron gas was not yet investigated within 
this approach. 

In the present work, the ground-state properties of a translation invariant A^-polaron 
system are theoretically studied in the framework of the variational path-integral method for 
identical particles, using a further development of the model introduced in Refs. 2^ 2^, 2^ . 
In Sec. II, the variational path-integral method and the chosen model system are described. 
In Sec. HI, we discuss the numerical results for the ground-state energy of a translation 
invariant iV-polaron system. Sec. IV is a summary of the obtained results with conclusions. 
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II. VARIATIONAL PATH-INTEGRAL METHOD FOR A iV-POLARON SYSTEM 



A. The many-polaron system 

In order to describe a many-polaron system, we start from the translation invariant A^- 
polaron Hamiltonian 

N 2 ^ N N 2 ( ^ \ 

i=l j=l l=\,+j oo I J k \i=l k / 

where m is the band mass, e is the electron charge, culo is the longitudinal optical (LO) 
phonon frequency, and Vfc are the amplitudes of the Frohlich electron-LO-phonon interaction 

= (^] ^'^ { ^ \ a = ^— f^!!^] (J- - 1.) (2) 

of course with the electron-phonon coupling constant a > 0, the high-frequency dielectric 
constant > and the static dielectric constant eg > 0, and consequently 

^>n('-^Y\^.V-2<(J^f.u, (3) 

Coo \ m J \nhJhoJ 

which is an important physical condition on the relative strength of the Coulomb interaction 
as compared to the electron-phonon coupling, as stressed in the earlier bipolaron work ji^l . 
In the expression (jSj), H* is the effective Hartree 

H* = ^, a*B = ^^ (4) 

where a% is the effective Bohr radius. The partition function of the system can be expressed 
as a path integral over all electron and phonon coordinates. The path integral over the 
phonon variables can be calculated analytically [i^. Feynman's phonon elimination tech- 
nique for this system is well known and leads to the partition function, which is a path 
integral over the electron coordinates only: 

Z=(t] ^f7r ]<fe'Vf (5) 

where r = {ri, • ■ ■ , riv} denotes the set of electron coordinates, and §T>y denotes the path 
integral over all the electron coordinates, integrated over equal initial and final points, i.e. 

//•f(/3)=r 
dv I e^Vf{T). 
Jr(0)=r 

Throughout this paper, imaginary time variables are used. The effective action for the 
A^-polaron system is retarded and given by 




N / r / \\ 2 ^ N N 9 



dr 



j=i i=i,^j \ J \ / 



9 w / ^ I \\ cosh hjji^o (^3 — \t — a\) 
2g«k.(r,(r)-n(a)) ^^M2A^ \ ^^dadr. 6 

sinh ^(3hioi^o 
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Note that the electrons are fermions. Therefore the path integral for the electrons with par- 
allel spin has to be interpreted as the required antisymmetric projection of the propagators 
for distinguishable particles. 

We below use units in which h = 1, m = 1, and cjlo = 1- The units of distance and 

1 /2 

energy are thus the effective polaron radius [h/ (mc^Lo)] and the LO-phonon energy hwi^Q. 



B. Variational principle 

For distinguishable particles, it is well known that the Jensen- Fey nman inequality [Ij, |2f 
provides a lower bound on the partition function Z (and consequently an upper bound on 
the free energy F) 

(7) 

for a system with real action S and a real trial action S'o-The many-body extension (Ref. [l9l |. 
p. 4476) of the Jensen- Feynman inequality, discussed in more detail in Ref. jl^], requires (of 
course) that the potentials are symmetric with respect to all particle permutations, and that 
the exact propagator as well as the model propagator are defined on the same state space. 
This means that both the exact and the model propagator are antisymmetric for fermions 
(symmetric for bosons) at any time. The path integrals in ((7j) thus have to be interpreted in 
terms of an antisymmetric state space. Within this interpretation we consider the following 
generalization of Feynman's trial action 




I) 



sinh 



with the variational frequency parameters v,w,uj. Because the coordinates of the fermions 
enter Eq. only through the differences rj (r) — r; (a), this trial action is translation 
invariant. 

Using the explicit forms of the exact © and the trial actions, the variational inequality 
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(jHI takes the form 



F(/3|iVt,iVj<Fo(/3|iV^,iV^) + — / ( 




.f3 / N 



w V 



AN(3 




sinh \(3w 
coshc^LO -\r - o-|) 



dadr 



dadr. 



sinh |/?t^LO 

(10) 

and it is clear that the minimization automatically implies > w'^. 

In the zero-temperature limit {jS oo), we arrive after some lengthy algebra at the 
following upper bound for the ground-state energy {N^, Ni) of a translation invariant 
A^-polaron system 

with 

E,,r (iVt, ^il^, ^) = ^ - + (A^i) + (iVi) 

+ Ecu (^t) + Ec\\ (iVi) + Ecu N{) 

+ K|| (iVt) + (iV^) + (A^T> A^i) > (11) 

where (A^) is the energy of N spin-polarized fermions confined to a parabolic poten- 
tial with the confinement frequency a;, Ec\\ [^Nui)) is the Coulomb energy of the electrons 
with parallel spins, Ec^i^^i ^i) is the Coulomb energy of the electrons with opposite 
spins, Ea\\ (A^T(i)) is ^1^^ electron-phonon energy of the electrons with parallel spins, and 
Eau i^h-^l) is electron-phonon energy of the electrons with opposite spins. The key 
steps in the derivation and the resulting analytical expressions for the terms of Eq. can 
be found in the Appendix. 



III. DISCUSSION OF RESULTS 

In the present section we summarize and discuss the main results of the numerical mini- 
mization of E^ar (^t' ^5 ^) with respect to the three variational parameters v, w, and 
u. The Frohlich constant a and the Coulomb parameter 

an = —F= = with - = — (12) 

^ 1-7] T] Ec ^ ^ 



-oo 



characterize the strength of the electron-phonon and of the Coulomb interaction, obeying 
the physical condition a > ao [see (jS))]. The optimal values of the variational parameters 
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v,w,and u are denoted fop,Wop, and cUop, respectively. The optimal value of the total spin 
was always determined by choosing the combination {N\^ Ni) for fixed N = + A?"! which 
corresponds to the lowest value (N) of the variational functional 

(N) = mm E^ar {N^,N - N^\vop, Wop, Uop) ■ (13) 

In Figs. 1 to 3, we present the ground-state energy per polaron (panel a), the addition 
energy (panel b), the optimal values of the variational parameters (panel c) and the total 
spin (panel function of the number of polarons. The addition energy is determined 

by the formula 

A (N) = E°{N+1)- 2E° (N) + E° {N - 1) . (14) 

In Fig. 1 we consider a highly polar system with a = ao = 7. The optimal value Uop (see 
panel c) for the confinement frequency uj is strictly positive (at least for N < 31). Therefore, 
the results of Fig. 1 are related to multipolaron states analogous to those investigated in 
Ref. ji^. This interpretation is confirmed by the fact that (see panel a of Fig. 1) the 
ground-state energy per polaron for = 2 is lower than that for = 1. For N > 2, 
the ground-state energy per polaron is an increasing function of A^, which means that the 
effective electron-phonon coupling weakens due to screening when the number of polarons 
increases. 

The addition energy (panel b of Fig. 1) oscillates, taking local maxima at even A^ and 
local minima at odd A^. This oscillating behavior reflects the trend of a stable multipolaron 
state to have the minimal possible spin. This trend is an analog of the pairing of electrons in 
a superconducting state. For even A^ (see panel d of Fig. 1) the total spin S is equal to zero. 
For odd A^, one electron remains non-paired and S = 1/2. Therefore, one intuitively expects 
that the states with S = are energetically favorable as compared to states with S = 1/2, 
and hence, A (A^) for odd A^ is lower than A (A^) for even A^.The plot of the addition energies 
in panel b of Fig. 1 confirms this expectation. Furthermore, pronounced maxima in A (A^) 
correspond to closed-shell systems with A^ = 2,8,20.... 

The optimal values of the variational parameters (panel c of Fig. 1) reveal a general trend 
to decrease as a function of A^, except the parameter v, which have a peak at A^ = 2. This 
peak, as well as the minimum of E^ (N) /N a.t N = 2, shows that the two-polaron state in the 
extremely strong-coupling regime is especially stable with respect to the other multipolaron 
states with N > 2. The dependence of the parameter u on N starts from N = 2, because 
the one-polaron variational functional does not depend on u. 

In Fig. 2, the ground state energy, the additional energy, the variational parameters, 
and the total spin for A^-polaron systems are plotted for a = 3, ao = 4.5, and 1] = 1/3. 
In this regime, the optimal value for the confinement frequency u is Uop = (panel c). 
Therefore, in this regime, as well as at weaker electron-phonon coupling strengths, A^ > 1 
polarons do not form a multipolaron state. The addition energy, as seen from panel b of 
Fig. 2, has no oscillations or peaks in the case when A^ polarons do not form a multipolaron 
state. It should be noted, that in the case when Uop = 0, we deal with a finite number A^ of 
polarons in an infinite volume. So, at Uop = the many-body effects, related to the electron- 
electron interaction and to the Fermi statistics, are vanishingly small. The dependence of 
the ground-state energy of the total spin of a many-polaron system is just one of these 
many-body effects. As a consequence, the ground-state energy within the present model at 
ujop = does not depend on the total spin. For this reason, there is no panel d in Fig. 2. 

Figs. 1 and 2 represent two mutually opposite cases (with u ^ and with ujop = 
for all A^ under consideration). Fig. 3 describes the case when the regime with Uop 7^ 
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(for < 16) changes to the regime with uoop = (for > 17). As seen from panel 
a of Fig. 3, the ground-state energy for A^ < 16 behaves similarly to that calculated 
for a = ao = 7 (panel a of Fig. 1), with the following distinction: for ao/a = 1.01 
(a = 7) it appears that E'^ (N) /N\j^^^ > E'^ (N) /N\j^^^, while for ao/a = 1 {a = 7), 
E'^ (N) /N\j^^2 < (^) /^lAr=i- As seen from panel c of Fig. 3, when an extra polaron 
is added to A^ = 16 polarons, the optimal value for u switches to zero, and therefore, the 
multipolaron state transforms to the ground state of A^ independent polarons. When A^ 
changes from A^ = 16 to A^ = 17, the ground-state energy per polaron slightly jumps down 
and is practically constant with further increasing A^. The transition from a multipolaron 
state to a state of A^ independent polarons is clearly revealed in the dependence of the 
addition energy on the number of polarons (panel b of Fig. 3). At A^ = 16, A (A^) has a 
pronounced minimum, which is a manifestation of the transition from a multipolaron ground 
state to a ground state of N independent polarons. 

The total spin, as seen from panel d of Fig. 3, takes its minimal possible value for A^ < 13, 
while for A^ > 14, the ground state is spin-polarized. So, the transition from the ground 
state with the minimal possible spin to the spin-polarized ground state with increasing A^ 
precedes the break-up of a multipolaron state. For A^ > 17, in the same way as in the 
case (a = 3, ao = 4.5) , the variational ground-state energy of an A^-polaron system does not 
depend on the total spin. 

In Fig. 4, the "phase diagrams" analogous to that of Ref. are plotted for an A^- 
polaron system in bulk with A^ = 2, 3, 5, and 10. The area where < a is the non-physical 
region. For a > ao, each sector between a curve corresponding to a well defined A^ and 
the line indicating ao = a shows the stability region where Uop 7^ 0, while the white area 
corresponds to the regime with Uop = 0. When comparing the stability region for A^ = 2 
from Fig. 4 with the bipolaron "phase diagram" of Ref. [2^ , the stability region in the 
present work starts from the value ac ~ 4.1 (instead of etc ~ 6.9 in Ref. |2|^). The width 
of the stability region within the present model is also larger than the width of the stability 
region within the model of Ref. j27[. Also, the absolute values of the ground-state energy 
of a two-polaron s;^tem given by the present model are smaller than those given by the 
approach of Ref. j27| . 

The difference between the numerical results of the present work and of Ref. [2^ is due 
to the following distinction between the used model systems. The model system of Ref. 
(27! consists of two electrons interacting with two fictitious particles and with each other 
through quadratic interactions. But the trial Hamiltonian given by Eq. (6) of Ref. [2^ is 
not symmetric with respect to the permutation of the electrons. It is only symmetric under 
the permutation of the pairs "electron + fictitious particle". As a consequence, this trial 
system is only applicable if the electrons are distinguishable, i.e. have opposite spin. In 



contrast to the model of Ref. [27[ , the model used in the present paper is described by the 
trial action (9), which is fully symmetric with respect to the permutations of the electrons, 
as is required to describe identical particles. Up to now we have been unable to construct 
such a model with two retardation sources. As a consequence, the trial model of Ref. [27] is 
superior to our model for describing 2 polarons because it has more varaiational parameters, 
but its applicability is limited to 2 polarons. The generalization of the model of Ref. 113] to 
A^ > 2 is currently under investigation. 

The "phase diagrams" for N > 2 demonstrate the existence of stable multipolaron states 
(see Ref. (3). As distinct from Ref. [i^, here the ground state of an A^-polaron system is 
investigated supposing that the electrons are fermions. As seen from these figures, for N > 2, 
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the stability region for a multipolaron state is narrower than the stabihty region for N = 2, 
and its width decreases with increasing A^. The critical value etc for the electron-phonon 
coupling constant increases with increasing A^. From this behavior we can deduce a general 
trend, which explains the behavior of the ground-state energy and related quantities as a 
function of shown in Fig. 3. For fixed values of a and rj, the width of the stability region 
for a multipolaron state is a decreasing function of the number of electrons. Therefore, 
for any (a,?]) there exists a critical number of electrons Nq {a,r]) such that a multipolaron 
state exists for N < Nq {a.rj) and does not exist for N > Nq (a,?]). For example, for the 
A^-polaron system described in Fig. 1, A^o is at least larger than 20. For the system shown 
in Fig. 2, A^o = 1? and for the system in Fig. 3, A^o = 16. If we add electrons to an A^-polaron 
system one by one, the multipolaron state breaks up when the number of electrons exceeds 
a critical value A'o 

In order to analyze the consequences of the Fermi statistics for the ground-state properties 
of an A^-polaron system, we compare the ground-state energies calculated with and without 
the Fermi statistics. In Table 1, the results are presented for the ground-state energy per 
particle in units of the one-polaron strong-coupling energy Ei, 



with a = 10, = for three cases: the many-body path-integral approach of the present 
work with fermion statistics (S^^^), the same approach for distinguishable particles (S^^), 
and the strong-coupling approach of Ref. [HI, which also does not take into account the 

Fermi statistics {S^'^'^^). 



Table 1. The polaron characteristic energy Sn calculated using different methods. 



N 


c(F) 


cW 


Md,sc) 
<-'N 


2 


-1.349 


-1.349 


-1.148 


3 


-1.308 


-1.415 


-1.241 


4 


-1.296 


-1.468 


-1.308 


5 


-1.279 


-1.508 


-1.361 


6 


-1.272 


-1.536 


-1.404 



As seen from Table 1, the ground-state energy per particle for A^ identical polarons Sj^ 
is higher than that for A^ distinguishable polarons Furthermore, increases whereas 
S^^^ decreases with an increasing number of polarons. Note however that S^^ < S^'^^'^ for the 
considered values of a and t], which means that the path-integral variational method provides 
better results for the A^-polaron ground-state energy than the strong-coupling approach [i^ 
(at least for a < 10). 

Another consequence of the Fermi statistics is the dependence of the polaron character- 
istics and of the total spin of an A^-polaron system on the parameters (a, ao,A^). In Fig. 
5, we present the ground-state energy per particle, the confinement frequency and the 
total spin 5* as a function of the coupling constant a for ao/a = 1.05 and for a different 
numbers of polarons. The ground-state energy turns out to be a continuous function of 
a, while cUopand S reveal jumps. For all considered numbers of polarons N > 2, there is 
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a region of a in which S takes its maximal value, while uJop 7^ 0. When lowering a, this 
spin-polarized state with parallel spins precedes the transition from the regime with lOop 7^ 
to the regime with ujop = (the break-up of a multipolaron state). For N = 2 (the case of a 
bipolaron), we see from Fig. 5 that the ground state has a total spin S = for all values of 
a, i. e., the ground state of a bipolaron is a singlet. This result is in agreement with earlier 
investigations on the large-bipolaron problem (see, e. g., |3l|). 

IV. CONCLUSIONS 

Using the extension of the Jensen- Feynman variational principle to the systems of iden- 
tical particles, we have derived a rigorous upper bound for the free energy of a translation 
invariant system of interacting polarons. In the zero-temperature limit, the variational 
free energy provides the variational functional for the ground-state energy of the A^-polaron 
system. The developed approach is valid for an arbitrary coupling strength. The resulting 
ground-state energy is obtained taking into account the Fermi statistics of electrons. 

For sufficiently high values of the electron-phonon coupling constant and of the ratio 
= Eo/Eooi the system of interacting polarons can form a stable multipolaron ground 
state. When this state is formed, the total spin of the system takes its minimal possible 
value. The larger the number of electrons, the narrower the stability region of a multipolaron 
state becomes. So, when adding electrons one by one to a stable multipolaron state, it breaks 
up for a definite number of electrons Nq, which depends on the coupling constant and on the 
ratio of the dielectric constants. This break-up is preceded by the change from a spin-mixed 
ground state with a minimal possible spin to a spin-polarized ground state with parallel 
spins. 

For a stable multipolaron state, the addition energy reveals peaks corresponding to closed 
shells. At N = No, the addition energy has a pronounced minimum. These features of the 
addition energy, as well as the total spin as a function of the number of electrons, might be 
resolved experimentally using, e.g., capacity and magnetization measurements. 
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APPENDIX A: MATHEMATICAL DETAILS 

1. Generalization of the Hellman- Feynman theorem 

For the averages of the quadratic terms in Eq. (jiop . we can derive a generalization of the 
Hellman-Feynman theorem for the case where we have a (trial) action but no Hamiltonian. 
Indeed, since Fq = —^InZo it follows that 

d _ 1 d _ 1 1 d _ 1 /dSo\ 
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for any parameter 7 in the trial action. Taking the derivative of 5*0 [eq. (jH))] with respect to 
uo and v then gives 



2iV /c/SoX 2N/3dF^ 



\J,(=i / \ / 







and therefore the variational inequality becomes 

+ 1/ (,4|r,(r)-r,M|)/- 

./3 / 



2/5 70 Jo V /o sinhl/?cuLo 

(A.2) 



2. Correlation and density functions 

In order to calculate the Coulomb and the electron-phonon energies [the terms in the 
second and third lines of Eq. (jA.2j) . respectively], we only need the pair correlation function 
qf and the two-point correlation function Cp for fermions which we define as 



gp {r,P\N^,N^) 



N{N 



1 ^ 



1 ^ 



(A.3) 
(A.4) 



where (. . .)q denotes a path-integral average with the action functional Sq. After a separation 
of the center-of-mass motion (see Ref. l25|), these correlation functions take the form 



^^(r,/3|Ar^,Arj = l_L_ ^ (5 (r - r, + r,))^ , 

/ w^r (/5 — r) 



(g,r, P\Nh Ni) = Cp {q,T, P\N^, N^) exp 



+ 



— w"^ sinh (|f r) sinh (if (/3 — r)) 



sinh 

sinh (|i-i^t) sinh (/? — r)) 



UJ 



sinh (leu/?) 



(A.5) 



(A.6) 



10 



with 

N 



CA^,r,P\N^,N,) = ^ (e-'i-.Me^'i-r.W)^, (A.7) 



iV2 

where {■ ■ ■)p denotes a path-integral average with the action functional 

r-/3 ^ 



i=i L 



dr 



-^'r] (r) 



dr (A.8) 



for = A^^i + A^^i independent fermions in a 3D parabolic potential with the confinement 
frequency uo. We shall use also the density function 

1 ^ 

The functions ()A.5j) . ()A.7jl . and ()A.9|1 were already derived before (see Refs. IHHIaS). 

Both the Coulomb energy and the electron-phonon energy in Eq. ()A.2|) are effectively 
Coulomb terms but with two important differences. Firstly, the standard Coulomb repulsion 
between the electrons is static, whereas the effective Coulomb attraction due to the polaron 
effect is retarded. A direct consequence of this difference is that the center of mass plays no 
role in the Coulomb repulsion, whereas it is essential in the retarded contribution. Secondly, 
the self-interaction has to be excluded from the Coulomb repulsion, whereas it contributes 
in the electron-phonon contribution. This is the main reason why we treat the Coulomb 
repulsion via the pair correlation function (in real space), and the retarded interaction with 
the two-point correlation function Cj? (k,r, /3|A^) (i.e. in wave number space). In principle 
we have the choice to handle both terms either in real space or in wave number space. 

Having the definitions of the pair correlation function and the two point correlation 
function in mind, we thus obtain for the free energy (performing the angular integrations at 
once) 



F{f3\N^,N^)<Fo{(3\NhNi)--uJ 



1 dF,{(3\N^,N{) lv^-w^dFo{(5\N^,N{, 



2 du 2 V dv 

POO 

+ 2txU / rN {N -l)gF{r,(3\N^,N{)dr 



TT 



/•p/^ r°° coshM-0 — r 

I I N'CF{q,T,mhNi)dq^^^^^dr. (A.IO) 



For the ground state energy (/5 oo) we thus find 

E° (AT^, N^) < E^ar {N^,N^\v, w, to) , 
where the variational functional is 



2 duj 2 V dv 

+ 2-nU I rN{N - 1) gp (r, f3 ^ oo|A^|, Ni) dr 



V2a 



TT 



/•oo /•oo 

/ / N^CFiq,T,(3 oo\Ni;,Ni)dqe-^dT. (A.ll) 
Jo Jo 
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Here, Eq is the ground state energy corresponding to the trial action, given by 

Eo {N^, N^) = ^{v-w-uj) + Ef (N^) + Ep (N^) , (A. 12) 

where Ep {N) is the ground state energy of fermions with parallel spins and with energy 
levels e„ = (n + |) cu: 

Ep (N) = Il{L + 2) (L + lfu; + {N- N^) (^L + cu, (A.13) 
L denotes the lowest partially filled or empty level, and 

Nl = ]:L{L + 1) (L + 2) 
o 

is the number of fermions in the fully occupied levels. Filling out Eq in E^ar we thus obtain 

3 ( v uj^^ 3 1 1 
E^ar {N^,N^\v,w,u) = - ^""""^^ --UJ + -Ep (ATj + -E^ (N^) 



+ 2nU I rN{N- 1) gp (r, /S oo|A^|, A^^) dr 
V2a 



TT 



POO PQC 

/ / N'^Cp{q,T,p oo\N^,N^)dqe-^dT, (A.14) 
Jo Jo 

where the factor 1/2 in front of Ep is a consequence of the subtraction Ep — ^i^^f, not sur- 
prising because of the virial theorem for harmonic oscillators. The term | is precisely 
the same as in Feynman's treatment of the polaron, but of course the values of v and w will 
be different if many-particle effects will be taken into account. 

We now split gp and Cp for a mixture of fermions with different spin projections into 
the parts corresponding to parallel and opposite spins. The case of electrons with spin 
up and A^j^ electrons with spin down can be found after some reflection in terms of the 
spin-polarized quantities: 

Nfip {q,P\N = N^ + N^) = N^np (q,/3|Ar^) + N^fip (q,/5|Arj , (A.15) 
N{N-l) gp (r, /?| AT^, N{) = {N^ - 1) gp (r, P\N^) + A^^ {N^ - 1) gp (r, P\N^) 

+ / e^''"^F(q,/5|A^T)^F(q,/?|Ari)rfq, (A.16) 

(27r) J 

N^Cp (q,r, P\N^, N^) = NfCp (q,r, P\N^) + NfCp (q,r, P\N^) 

+ 2N^N^np {q,md np (q,/5|Arj . (A.17) 

3. Coulomb and electron-phonon energies 

Using Eqs. ()A.15|) to (jA.17|l in the variational functional (jA.llll . we arrive at the expres- 
sion with three Coulomb terms and three electron-phonon terms as follows: 

E,,r (A^t, AT^, u) = \^^^^^ + {N{) + h.p (ATJ 

+ Ec\\ {N^) + Ec\\ {N{) + Ecu N{) 

+ K|| {N^) + {N^) + N^) (A.18) 
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with the contributions 

Ec\\{N) = 2txN {N -l)U j rgFir,P ^ oo\N)dr, (A.19) 

Jo 



(A.20) 



Ea\\ (N) 



a 



71 




2 / 2 2 2 / \ 1 —LJT \ 



^0 



oo /"OO 



(ire 



Ci? (q,r, (5 oo|A^) dqdr, 
(A.21) 



JO JO 

X hp (q,/? ^ oo|A''-f) riiT (q,/? ^ oolA?"!) . 



(A.22) 



The integrations over q and over r in Eqs. ()A.19|) to (|A.22p are performed using the exphcit 
form of the density and correlation functions (see Eq. 25]) 



1 ^ 

hp (q,/3 ^ oojAT) = - ^ (g) A (A;|/3 ^ oo, AT) 



(A.23) 



fc=0 

L oo 



(q,r, /5 - oo| AT) = ^ ^ M^y (q) e^'-'>^ [A (A;|/3 ^ oo, AT) - (A;, A;'|/3 ^ oo, AT)] 

(A.24) 



k=0 k'=L 



L oo 



1 r ^ ^ 



N{N-1) {27ry 
X [/i {k\P ^oo,N)- h {k, k'\P ^ oo, N)] 



with the matrix elements 

(g) = exp (|- 



(A.25) 



(A.26) 



Mkk' (g) = e" 



2cj 



/ ^2 \ fc>-fc< fc< 



(j + 1) 



(fc<-j)! _ 



(fc>-fc<) / g 



A;< = min {k, k') , 
A;> = max {k, k') J ' 

(A.27) 



where L^""* (x) are the Laguerre polynomials, and with one-particle and two-particle distri- 
bution functions 

1, k < L, 

fi{k\l3 ^ oo,N) = { 0, k>L, (A.28) 

^^=^ k = L, 



h{k,k'\(3^oo,N) 



' /i {k\i5 ^ oo, AT) /i ^ oo, AT) , A; ^ A;', 

1, k = k' <L, 

0, k = k' > L, 

N-Nl N-Nl-1 U ^ U ^ t 

Nl+i-NlNl+i-Nl-V a. a. i^. 



(A.29) 
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Using Eqs. ()A.23|) to ()A.29|1 . after performing integrations we arrive at the following 
formulae for the Coulomb and electron-phonon energies ()A.19|) to ()A.22|) . 
(i) The Coulomb energy for opposite spins is given by the expression 



k=0 1=0 ^ 



k + l-^ 



X 



Li + 2 



+ 



Ni-NL 



Li + 2 



1 + 3 J Nl,+i -Nl,\1 + 2 
(ii) The Coulomb energy for parallel spins is 

I L L 

Ec\\ [N) 



(A.30) 



k=0 1=0 ^ 



I 



+ 
+ 



L + 2\/L + 2\ N-Nl /L + 2\ /L + 2 
k + sjyi + s) NL+i-NL\k + 2j\l + 3 
N -Nl fL + 2\fL + 2 



Nl+i- Nl\1 + 2 J \k + 3 
N-Nl N-Nl-1 fL + 2\fL + 2 



Nl+1 - Nl Nl+1 

L L 



* k=0 k'=0 



k< I / ■ 1 



A;> - A;< + / + j - i\ fky + 2\ [2 {I + ky - k 



1=0 j=0 

(iii) The electron-phonon energy for opposite spins is 

Lt L 



l-J 



-2a 



dre^Y^Y} ^^^^ ^ ^ 



Lt + 2 



+ 



,=o«=o {2ujP{t) + \\ 
- Nl, (L^ + 2 



k+l- 



k + 3j ■ Nl,+i 



NL.\k + 2 



Li + 2 
1 + 3 



+ 



N^-Nl, (L^ + 2 



with the function 

P(r) 



w 

^r + 



2{N^ + N{) 



Nl,\1 + 2 



(1 - e"^") 



(A.31) 



(A.32) 



(iv) Finally, the electron-phonon energy for parallel spins takes the form 

V2a 



{N) = j dre-^E^i^ {N,t) 



(A.33) 
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where the time-dependent function E^w {N, r) is a sum of three terms: 

E., iN,r) = E^, iN,r) + (iV,r) + _ _ _ (^'-) ■ 

(A.34) 

The terms ^Jy^ (iV, r) can be written down in two equivalent alternative forms. The first 
form is relevant for the numerical calculation in the region of small and intermediate values 
of (a;r), 




L-l L-j-l n / I \n+m (' L+2 \ /20+n)Wrn-iwj+n+m-i\ 

/ jujT _|_ -ju}T\ V / Vj+n+3/ V n-m )\m>\ j+n ) 

h ' [2a,P(r) + ir"+"+i 

(A)GS)(2!l('""*) 



L-l L-l / I \j+n /L+2\ (L+2\ (j+n-^\ 

yy^_^ [m!ll^±^l\^J^_^ \ (A.35) 



^ [2(cosha;r-l)]^ (^+3e7) 



'^"^ ^ V 2 [2(^P(r) + l-e— 

L L—j n 



i'-i n / -i \n+m / L+2 \ p(i+n)Wm-|\ /'j+n+m-|\ 
^ ' \j+n+2/ V n-m / V m / \ j+n ) 



„.o™.o [2a.P (r) + 

' (-ir "(ta(asp ?)rr^) 



+2 E E ,„ .'^ii. " 1 ■ 



all ^ ' ^ A/ o I /./.I ^ [2o;P (r) + 1^^^+^ 




j=0 n=0 



[2a;P(T) + lf+"+^ 
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The second form is relevant for intermediate and large values of (out), 

I / oo L-1 n /_i \n+m /'L+j+2W2(«+j)Wm-|wi+n+m-|\ 

^(0) ^ ^ i ^-jujT l Kn+j+sJ \ n-m A m A j+n ) 



L-l L-j-1 n 



> \n+j+V \ n-m )\m)\ j+n ) 



^j+n+m+\ 



_ g^i'^T V \n+i+3/ V n-m / V m / V j+r; 

j=i 71=0 m=o [2c<jP (r) + l]"* 

I y^ V I \n+-i) lm+3/ I n A m / \ (A 38) 

^£-0 [2^P(r) + ir-+^ 



oo L n / -r Vfi + 'fi /L^;/>2\ /2(r/+,/)Wm-iwj+7i+7fi-i\ 

L-j n yn+m / L+2 W2(n+i)Wm-|\ p+n+m-i\ 

^ ' VTi+7+2/ V n-m A m A i+n / 



Z—^ Z_-^ ro n / \ , i]J+n+m+-5 



•i+j+2/ V n—m /V m /V j'+n 

7=1 n=0 m=0 



L L-l (_-\\n+m rL+2\rL+2\rn-^\(n+m-^\' 
_|_2 y^ y^ V ' Vn+2Am+3A n )\ m ) \ ggN 

[2a;P(r) + ir-+^ 
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Fig. 1 



Fig. 1. (Color online) The ground-state energy per polaron (a), the addition energy (6), 
the optimal values of variational parameters (c) and the total spin {d) as a function of 
for a translation invariant A^-polaron system with a = 7, ao = a. 
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Fig. 2. (Color online) The ground-state energy per polaron (a), the addition energy (b), 
and the optimal values of variational parameters (c) as a function of for a translation 
invariant A^-polaron system with a = 3, ao = 1.5a. 
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Fig. 3 



Fig. 3. (Color online) The ground-state energy per polaron (a), the addition energy (b), 
the optimal values of variational parameters (c) and the total spin {d) as a function of 
for a translation invariant A^-polaron system with a = 7, ao = l.Ola. 
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Fig. 4 



Fig. 4. (Color online) The "phase diagrams" of a translation invariant A^-polaron system. 
The grey area is the non-physical region, for which a > a^. The stability region for each 
number of electrons is determined by the equation ac < a < ao- 
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S. N. Klimin, F. Brosens, and J. T. Devreese 



Fig. 5. (Color online) The ground-state energy per particle (a), the optimal value uJop of 
the confinement frequency (6), and the total spin (c) of a translation invariant A^-polaron 
system as a function of the coupling strength a for ao/a = 0.5. The vertical dashed lines 
in the panel c indicate the critical values ac separating the regimes of a > ac, where the 
multipolaron ground state with Uop 7^ exists, and a < ac, where Uop = 0. 
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